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Discontinuous Galerkin methods for hyperbolic and 
advection-dominated problems on surfaces 
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We extend the discontinuous Galerkin (DG) framework to the analysis of first-order hyperbolic and 
advection-dominated problems posed on implicitely defined surfaces. The focus will be on the hyper¬ 
bolic part, which is discretised using a “discrete surface” generalisation of the jump-stabilised upwind 
flux considered in |Brezzi et al. |f2004) . A key issue arising in the analysis (which does not appear in the 
planar setting) is the treatment of the discrete velocity field, choices of which play an important role in the 
stability of the scheme. We then prove optimal error estimates in an appropriate norm given a number of 
assumptions on the discrete velocity field, which are then investigated and discussed in more detail. The 
theoretical results are verified numerically for a number of test problems exhibiting advection-dominated 
behaviour. 

Keywords : discontinuous galerkin; upwind; surface partial differential equations; hyperbolic partial dif¬ 
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1. Introduction 


Partial differential equations (PDEs) on manifolds have become an active area of research in recent 
years due to the fact that, in many applications, mathematical models have to be formulated not on 
a flat Euclidean domain but on a curved surface. For example, they arise naturally in fluid dynam¬ 
ics (e.g., surface active agents on the interface between two fluids,[James & Lowengrub ( 20043 ) and 
material science (e.g., diffusion of species along grain boundaries, Deckelnick et al. ( 2001] >) but have 
also emerged in other areas as image processing and cell biology (e.g., cell motility involving processes 
on the cell membrane, Neilson et aTj (j20TTji or phase separation on biomembranes, jElliott & Stinner 

( |20T0l >). 

Finite element methods (FEMs) for elliptic problems and their error analysis have been successfully 
applied to problems on surfaces via the intrinsic approach in |Dziuk| ( |1988} >. This approach has subse¬ 
quently been extended to parabolic problems |Dziuk & Elliott (2007b i as well as evolving surfaces Dziuk 
& Elliott (2007a). The literature on the application of FEM to various surface PDEs is now quite exten¬ 
sive, a review of which can be found in Dziuk & Elliott ( 201 3[ >. High order error estimates, which require 
high order surface approximations, have been derived in Demlow (2009) 1 for the Laplace-Beltrami op¬ 
erator. However, there are a number of situations where conforming FEMs may not be the appropriate 
numerical method, for instance, problems which lead to steep gradients or even discontinuities in the 
solution. Such issues can arise for problems posed on surfaces, as in Sokolov et al. < 2012) ) where the 
authors analyse a model for bacteria/cell aggregation. Without an appropriate stabilisation mechanism 
artificially added to the surface FEMs scheme, the solution can exhibit a spurious oscillatory behaviour 
which, in the context of the above problem, leads to negative densities of on-surface living cells. 

Given the well-known in-built stabilisation mechanisms discontinuous Galerkin methods possess 
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for dealing with hyperbolic/advection dominated problems and solution blow-ups, it is natural to extend 
the DG framework for PDEs posed on surfaces. DG methods have first been extended to surfaces in 
Dedner et al. ( 2013) >, where an interior penalty (IP) method for a linear second-order elliptic problem 
was introduced and optimal a priori error estimates in the L 2 and energy norms for piecewise linear 
ansatz functions and surface approximations were derived. These results were then generalised to both 


a more general class of surface DG methods and higher order surface approximations in Antonietti et al. 
( |20T4| . 

In this paper, we will extend the analysis of the surface DG method to the model problem 


—eAru + 'Vr ■ (wu) + cu = / onT, 


( 1 . 1 ) 


where e is assumed to be small or even equal to zero. For the elliptic part (when present) we can make 
use of the surface DG framework described in Antoniett i et q (.]( |2014| >. Our main focus will thus be on 
the hyperbolic part. The corresponding model problem takes the form 

Vj- • (wu) + u = f onT, 

where F is a compact smooth oriented surface in R 3 , c > 0 and w is a velocity-field which is purely 
tangential to the surface F. This advection problem is discretised using a “discrete surface” generali¬ 
sation of the jump-stabilised upwind flux considered in Brezzi et al. ( 2004j >. A number of challenging 
issues which do not appear in the planar setting arise when attempting to prove stability of the numerical 
scheme, related to the treatment of the velocity field on the discrete surface. We derive optimal a priori 
error estimates for this scheme given a number of assumptions on the discrete velocity field. We then jus¬ 
tify these assumptions by choosing the discrete velocity field to be a Raviart-Thomas-type interpolant of 
the velocity field. Numerical results are then presented for test problems exhibiting advection-dominated 
behaviour, suggesting that our surface DG method is stable and free of spurious oscillations. 


2. Notation and setting 

2.1 Continuous surface F 

Let r be a compact smooth and oriented surface in R 3 given by the zero level-set of a signed distance 
function |c/(x)| = dist(x,T) defined in an open subset U of R 3 . For simplicity we assume that dr = 0 
and that d < 0 in the interior of F and d > 0 in the exterior. The orientation of F is set by taking the 
normal v of F to be pointing in the direction of increasing d whence 


v(^) = v^Uer. 

With a slight abuse of notation we also denote the projection to F by t;, i.e. ^ : U —> F is given by 

B,(x) =x — d(x)v(x) where v(x) := v(t;(x)). (2.1) 

It is worth noting that such a projection is (locally) unique provided that the width Su > 0 of U satisfies 


8u < 


m ax||'q|lL“(r) 

1 = 1,2 


where Kj denotes the zth principle curvature of the Weingarten map H, given by 

H(x) := V 2 d(x). 


( 2 . 2 ) 
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Later on, we will consider a triangulated surface f/, C U approximating F such that there is a one-to-one 
relation between points x G f/, and % G F so that, in particular, the above relation (2.1 1 can be inverted. 
Throughout this paper, we denote by 


p($):=i -v(«)®v(?Uer, 

the projection onto the tangent space T^F on F at a point | G F. Here 0 denotes the usual tensor 
product. 

Definition 2.1 For any function q defined on an open subset of U containing F we can define its 
tangential gradient on F by 

V r rj := Vt 7 — (V 77 • v) v = PV 77 
and then the Laplace-Beltrami operator on F by 


A r ri := V r ■ (V r rj). 


Definition 2.2 We define the surface Sobolev spaces 

H m (r):={u€L 2 (r) : D a u G L 2 (T) V|a| 0 m}, mG NU{0}, 
with corresponding Sobolev seminorm and norm respectively given by 

V 2 / 

' ( m 

> ll“l!// m (r) : = M/F(n 

\/t=0 

We refer to |Wloka] ( j1987] ) for a proper discussion of Sobolev spaces on manifolds. 

Throughout this paper, we write x < y to signify x < Cy, where C is a generic positive constant whose 
value, possibly different at any occurrence, does not depend on the meshsize. Moreover, we use x ~ y 
to state the equivalence between x and y, i.e., C\y 0x0 Coy, for C 1 , C 2 independent of the meshsize. 


'I H m (r) I Yj W dCCu W l 2 (r) 



2.2 Discrete surface F/, 

The smooth surface F is approximated by a polyhedral surface F/, C U composed of planar triangles. 
Let ■% be the associated regular, conforming triangulation of f/, i.e. 

J* = U K h . 

K h etX h 


The vertices are taken to sit on F so that f/, is its linear interpolation. We assume that the projection map 
4 defined in (2.1 1 is a bijection when restticted to f/,, thus avoiding multiple coverings of F by I/,. Let 
< 0 /, denote the set of all codimension one intersections of elements . K h G 5;, (i.e., the edges). We 
define the conormal nt on such an intersection e/, G S\) r of elements K, and K7 by demanding that 


• nt is a unit vector, 

• nt is tangential to (the planar triangle) Kt, 

• in each point x G ei, we have that nt ■ (y — x) 0 0 for all y G . 

Analogously one can define the conormal n h on e/, by exchanging K h with Kf. It is important to note 
that 
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in general, and in contrast to the planar setting. Finally, we will denote by V/, the outward unit normal 
to !/, and define for each K b G ■% the discrete projection P/, onto the tangential space of I/, by 

P h ( x ): = I - V/, (x) <g> V/, (x), x G r h , 


so that, for v/, defined on f/ ; , 

Vjjv/, = P/iVv/,. 

Let fT C R 2 be the (flat) reference element and let F^ h : K K b C R 3 for K b G -f ,. We define the DG 
space associated to f/, by 

V h = {v h G L 2 (T b ) : v h \ Kh = X °F Kh l , X € P 1 (K) VK h G %}. 

For Vh G Vh we adopt the convention that vt is the trace of v/, on e/, = K h D K ti taken within the interior 
of Ar, respectively. 


2.3 Relating 17, to L 

Definition 2.3 For any function w defined on 1/, we define the surface lift onto F by 

V l {Z):=ri(x{Z)),ZGr, 


where by ( 2 . 11 and the non-overlapping of the triangular elements, x{£f) is defined as the unique solution 
of 

x = B, +d(x)v(t;). 


Extending constantly along the lines s i—>■ ^ +sv(q) we obtain a function defined on U. By (2.1 1 , 
for every K b G there is a unique curved triangle K[ := ^ (Kh) C F. Note that we assumed |(jr) 
is a bijection so multiple coverings are in fact not permitted. We now define the regular, conforming 
triangulation G7} of F such that 

r= U K l 


The triangulation of F is thus induced by the triangulation .77, of 1], via the surface lift. The appro¬ 
priate function space for surface lifted functions is given by 


v h ; = M e L 2 (F) : v{,(|) = v,,(*(£)) with some v h GV h }. 


We also denote by 1 ] 1 G Vh the inverse surface lift of some function 77 G Vj r satisfying (rj l ) ! = 77 . 
Finally, by applying the chain rule for differentiation on ( |2.1[ ), one can show that for x G 1], and v b 
defined on 1 ],, we have that 


Vr A Vft(x) = Ph(x)(I- dH)(x)P(x)V r Vh(^ (x)). (2.3) 

Finally, for x G 1/,, we denote the local area and local edge deformations when transforming T], to F by 
respectively S/,(x) and S e/i (x) i.e. 

8h(x) dA h (x) = dA(^(x)), 8 ei fx) ds h (x) = ds(|(x)). 

We finally state and prove some geometric estimates relating f/, to F. 
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Lemma 2.1 Let T he a compact smooth and oriented surface in R 3 and let I), be its linear interpolation. 
Then, omitting the surface lift symbols, we have that 


IMIIi“(ii) ~ ^ 2 > (2-4) 

U-Sh\\ L ~(r h )<h 2 , (2.5) 

l|v-v*|| L ~(r h )<h, (2.6) 

P-4JI L-(* k )<h 2 , (2.7) 

\\n + l--Yn + h / -\\ L ^ h) <h 2 ( 2 . 8 ) 

||P-RP*|| L ~(r h )<h (2.9) 

||P-PP /7 P||^ (J7!) < /z 2 . (2.10) 


Proof. The hist three estimates are proven in Dziuk ( 1988[ >. The next two estimates are proven in 


Antonietti et al. 


(2014). For the last estimate, it is sufficient to show that (P — PP/,P)jr for x G R- 1 scales 
appropriately. Setting x = Px (which is tangential to F) and noting that P/,x = x — (x ■ V/,)V/,, we have 


that 


(P -PP,,P)x = x-PP h x = x-P(x- (x-Vh)Vh) =X— {x-(x-V)V-(x-Vh){Vh - (v/, • v)v) 
= (x-V h )v h - (x-v h )(v h • v)v = (x-v h )(v h - v) + (x-v h )v (1 - (v h ■ v)) 

= (jc-(v* - v))(v /; - v) + (x-(v h - v))v (1 - (v* ■ v)) < \v h - v| 2 |x| + \\v h - v| 2 |x| 

</l 2 |x|</l 2 |x| 


where we have used the equality 1 — (v/, • v) = \ \ V/, — v| 2 and the geometric estimate (2.6 1 . The proof 
of the second to last estimate follows similar arguments. □ 

We complete this section by debiting an L? type projection operator for function in // 2 (Fi,): 


Lemma 2.2 Let r\ G H 2 {rf) and denote by TI hk r) the Lr projection of 77 onto 14 for k = 1 and piecewise 
constant functions on Ij, for k = 0. Furthermore, for notational simplicity, we dehne n /,77 := IT/ tl T]. 
Then, for sufficiently small h, we have that 


II 1 ? ^hk'n\\LP{K h ) ~ ^ lk+ ^W r \Ww k+ kP{K h )i 1 ^ P ^ °°j ^ — 0,1, 

\\P-n h ri\\ L 2 (dKh) <h^ 2 \\ri\\ H 2 {Kh) 

for each K k G ■%. 


Proof. The proof of both estimates follow from applying standard arguments on each K/, G ■% (which 
hold since each triangle K/, is planar). □ 


3. Problem formulation, discretisation and properties 

3.1 Model problem 

We will split our model problem o> into two parts: an elliptic part and a first-order hyperbolic operator 
which, when written in weak form over H 1 (F), are respectively given by 

J eV r n-V r vdA 
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and 


L 


—wu ■ Vp v + cuv dA 


where the velocity field w can be assumed to be purely tangential to the surface, i.e. w • v = 0 every¬ 
where, since any normal contribution would vanish when multiplied with Vfit. We will also assume, 
for simplicity, that the velocity field is divergence-free which, together with w ■ v = 0, implies that 
V f - w = 0. Finally the mass term is multiplied by a bounded function c > 0. The weak problem then 
reads: find u £ H l (F) such that 



wm) • Vj-v T cuv dA = 



Vveff^r). 


Existence and uniqueness of a solution u € H 2 (r) follows from standard arguments. 


(3.1) 


3.2 Discretisation of the hyperbolic operator 

Before we define a DG discretisation, we introduce the following discrete surface trace operators: 
Definition 3.1 For q € nK he ^- h Lr(dKi l ), {< 7 } and [ q] are given by 

{<?} : = ^(<7 + +0 ; W\ :=< 7 + -q~ on e/! e 4- 

For (j),h e [n Khe3 r h L 2 {dK h )} 3 , {</>;«} and [<j>;h} are given by 

{(j);h} := -((f> + ■ h + — (f>- ■ h~), [<j)',h] := 0 + -h + + -h~ on e/, € 4- 

Now we can define the discrete bilinear form for the advection operator: 

&h(u h ,Vh):= Y / ~ w hUh~ Vr h Vh + (c + Yi,)uhn dA h + Y / w*4['7i] ds h + ( tf Wh {u h ,v h ) (3.2) 
K h e3r h jKh e h e£ h Je h 

Wh is a discrete velocity field which will be explicitly related to w further down. We define the surface 
upwind flux w^uj, by 

Whui, = {wh\n h }{u h } + Z eh [u h ] 

where q f , h := \ | { vv;,;/!/,} |. For reasons that will be clear later on, the discrete mass perturbation coeffi¬ 
cient Yh is given by 


y h =max{-Vr h -w h ,--(c + Vr h -w h )} (3.3) 

Finally an additional term is added to take possible discontinuities in the discrete velocity field into 
account: 


^w h {uh,Vh)'■= Y [ l[ w D n h}{ui,Vh} 

e h eS h J ^ 1 


(3.4) 
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We make the following assumptions on the approximation properties of the discrete velocity field: 
if w £ [W 2 ’“(A^)p, then for each K/, £ we have 

\\PhW- l -w h \\ L ~ {Kk) <hP', (3.5) 

II' w h\\u°( Kk ) ~ h P2 (3.6) 

where p\,pi > 0. Note that assumption ( |3.6[ > implies that 

7 /, = —V^w/, for h small enough. (3.7) 


Remark 3.1 Note that 

_ [ KS w t' n h- w h_- n h)( v t- v h) if{ W h ;«/,}• n + >0; 

WhUh[vh\ = < u h (w,} • nX - w h ■ n h )(vj - v h ) if {w h \n h } ■ n j < 0; 

[ 0 if w h ■ nl = 0, 

and it can thus be seen that this flux works in exactly the same way as the classical (planar) upwind flux 
in the case that • nX — wf ■ nf = w ■ nf. 

We will discuss options for defining a discrete velocity field after proving stability and error esti¬ 
mates for the discretization. 


3.3 Stability of the advection operator 

In the following we shall prove stability of in the norm 


Wh II' Wl 2 (r h ) + I 


e h es h 


L 2 fh) 


We first require a useful formula which holds for functions in 

H 1 {3h) := {v\ Kh G H\K h ) : MK h £ £T h }. 

Lemma 3.1 Let 0 £ and yt £ Then we have that 

Y / Vf-nKb ds h = Y [<hn h \{Y} + {<l>',nh}[Vf}ds h . 

K h eSr h JdK h e h eS h Je h 

Proof. The result follows straightforwardly by noting that 

Y / ■ n K h ds h = Y / [^;«A]dsh= Y / [^;»/,]{W + {^;«/,}[v / l ds h- 

K h et? h JdK h e h eS h Je h e h eS h Je h 


Lemma 3.2 The surface DG bilinear form is stable, i.e., 


(3.8) 


□ 


for every £ V/,, provided that the discrete velocity field w/, satisfies (3.6 1 . 
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Proof. We proceed along the lines of Brezzi et al. <20041 by testing (3.2 1 with 17 , = w/, and integrating 
by parts on each K/, £ applying Lemma 3.1 By doing so, we obtain 


@ h {uh,u h )= Yj [ ( c + 'Yh + l v r h -w h ) u 2 h dA h + Y [ i w h'- n h }{«/,}[«/.] + &* IM | :2 dsh 

K h e3r h jK h \ 1 / e h eS h Je h 

^ {wh;n h }[u 2 h \ + [w h ;n h ]{u 2 h } ds h +- Y / {^h\n h ]{u 2 h } 

A e h e<g h Je h A e h eS h Je h 

= Y [ \ c +'Yh + l v r h -w h )u} l dA h 

K h e.% JK h \ Z / 

+ Y / {' W I>’ n h}{uii}[uh] + %e h I [uh] I 1 - \ { W h ; n h } [it*] 

e h eS h Je h A 

= Y [ ( c + Yh + lvr h -Wh) 1*1 d A h 

K h e3r h JKh \ A / 

+ Y* {Whinhjluhjluhl+Zejluh]] 2 - {w h -,n h }[u h ]{u h } 

e h eS h Je h 

= Y [ ( c + Yh + lvr h -Wh)u 2 h dA h + Y [ Ze h \[uh]\ 2 >\\uh\\ 2 H- 

K h e.9 h JK >’ V z / e h eS h Je h 


In the final estimate we used that c + //, + \ • w;, Js 

□ 


jc > 0 due to our definition of y/, given in (3.3 i. 


3.4 Discretisation of the advection-diffusion operator 
We consider a bilinear form 


Mi{uh,Vh) := %{uh,Vh) + &h{uh,Vh) 


(3.9) 


to discretize our problem O'- Here 3>h{uh-yh) is a discretization of —eAu on L/,. A number of 
formulations are given in Section 3.2 of Antonietti et af. il 20141. One such example is the symmetric 
surface interior penalty (IP) method: 


@h(uh,Vh)= Y / E ^h u h dA h + Y / 

K h e3r h jKh e h e£ h Je h 


ft{’i, [r/i] [R/i] ds h 
Y / ([ u h\{£Vr h Vh',n h } + [vh\{£'Vr h Uh-,n h }) ds h 

. Jeu 


(3.10) 


e h^h 


where [i eii := 5^ with a > 0 being a (penalty) parameter at our disposition. This method was first 
considered in Dedner et al. ( |2013) l for the elliptic problem including a mass term of the form Eitv. The 
method was shown to be stable for functions in V/, in the norm given by 


II I' 


Y I 


II 


I ||v& 


e h^h 


2 

L 2 {e h Y 


(3.11) 
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provided that a is large enough. With this in mind and taking Lemma 3.2 into account, we can easily 
prove the stability of in the following norm: 


112 .- 


II l 


— II' 11 z, 2 (J7, ) + e E 


K h e^h 




(3.12) 


3.5 Bilinear form on r 

We end this section by defining a bilinear form on F induced by SSh m . 



—wv i • Vj-n ’2 +CV 1 V 2 dA + 



({w;n}{vt} + ^[vi]) [v 2 ] ds. 


(3.13) 


where B, e i : = 8 eh 1 j | w ■ n | with n = n + or n =n is a surface conormals to e l h . Note that since w is assumed 

to be continuous on F and n + = —n we have \w;n\ = 0 and thus the addition term F lV/j vanishes on 
F and {w;n} =w -n + = — w ■ n~. Stability of &(•, •) in x follows from similar arguments as the 
proof of Lemma [F2] The norm associated with the bilinear form is given by 


II 


lilt 


IHlF(r)+ E 

e 'h eS i 



(3.14) 


Note the slight abuse of notation since here and in the following we use the same symbols to denote 
norms on F/, and on F. It will always be clear from the argument which norm we are refering too. 
Finally, we define 

^(vt,v 2 )= V / eV r vi • V r v 2 dA- V / lvi]{eV r v 2 ;n} + [v 2 ]{eV r vi;n} ds 

J “ » JK{ Jei 

+ E [ e Pe'l v lINds. (3.15) 


with Pei 8 e h l fie h - Note that is stable in Vj t x Vj t with respect to the DG norm given by 


:= e 


k h^h 




Ljlv^tw 


4 


2 

l 2 ( 4 ) 


(3.16) 


Furthermore it is bounded in // 2 (F) + x K- See 
Note that u satisfies 


Dedner et al. 


(2013) for further details. 


) := @(u,v) +&(u,v) = y I fv dA Vv G H 2 (T) + V/,. 
„,“L ,Jk! 




(3.17) 
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4. Convergence 


A DG discretisation of (3.11 is given as follows: find m/, € V/, such that 

.<dh(ph: Ui) / flA’h tlAh Vv/j € Vh 

K h eSr h jK h 

We will consider the error of the discretization in the norm 


lllr- 


II?, 


(4.1) 


on r. Before we state the main result, we make a note of key estimates relating norms defined on F to 
those on I/,: 

Lemma 4.1 Let v/, £ 14, then we have the following equivalence result: 


Furthermore, for u £H 2 (r), we have that 


,-h 


II tf 2 (F) 


< 


^lltFjr)- 


Proof. The proof of the first estimate follows similar arguments to that of Lemma 3.3 in Dedner et al. 
( j2013) >. The second estimate is given in (2.17) in |Demlow| ( |2009| . □ 

Theorem 4.1 Let u £ H 2 (r) and «/, £ V), denote the solutions to (3.1 1 and (3.21, respectively. Then, 
under assumptions ( |3.5| > and ( |3.6[ i on the discrete velocity field w;,, we have for h small enough: 

111“III ~ (lie 1 / 2 + h 3 / 2 +h Pl + s w (h + h Pl ~ 1 )^j ||w|| ff 2( r ) 


The constant s w is equal to 0 if | vi'/,= 0 and equal to 1 otherwise. 


We first give an outline of the proof: We split the error into two parts rj = u 1 — fli,u 1 and = 


,-i 


n h u 1 — u/, in the spirit of 


Brezzi et al. 


(2004): 


<-^lll<lll“-'-“/ ! |||<INIIH-|||zlll 


(4.2) 


The first term is a projection error and can be bounded using Lemma 2.2 and noting that |||r) Ill^lll 1 ? / Ill 
Lemma 4.2 Let rj be given as before. The we have that 

lll r llll M eZ +^ z )IMI H 2 {r) 

for h small enough. 

Using the stability of and £$/,(•, •) in V), x V/„ the second term in (|4.2[) can be estimated by 


lzllr^-*4(“ -Uh,x)-^h(ri,x), 


(4.3) 


Note that since we do not directly have Galerkin orthogonality, the first term on the right-hand side of 
( |4.3| > is not zero. We will discuss this term at the end of this section. We can deal with the second term 
in (|4.3|> using the following lemma: 
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Lemma 4.3 With 77 and % defined as before, we have that 

^h(ria)<H£' /2 + h 1/2 )\\u\\ H 2 ir) \\\x\\\ (4.4) 


for h small enough. 


Proof. The proof is a direct extension of the corresponding result in 
surfaces. We will first show that 


Brezzi et al. (2004ji to triangulated 


^h(r],x)<h 3/2 \\u\\ H 2 {r) 

Since Vj £ [V /,] 3 and P/, is constant on each K/, £ ■%, we have by definition that 

/ (n h oWh ■ V p h X ) V dA h = / (P h II h0 Wh ■ Vr h x) V dA h = 0. 

JK h JK h 

Using this, together with (3.7 1 , (3.6 1 , an inverse inequality on f/, and Lemma 2.2 the element integral 
term of 9%h with vi = 77 and V 2 = X becomes, for h small enough, 

E / -T? (wh-Vr h x) +( C + Yh)xv dA h = E / {( n hO' v h ~ w h ) ■ Vr h x) + 0 - v r h ■ w h )xv dA h 

uG.9J K h K[e3r h jKh 


K h eSr h 


~ E \\ n hOWh-Wh\\L°°(K h )\X\H l (K k ) + \\X\\L 2 (K h ) ) 1177 11 Z, 2 (AT*) 

<( E 1 \\X\\L 2 {K h ) + IIxHl 2 (K,,) ) N \\l 2 (Ki,) ~ ll M ll// 2 (r)IUIlL 2 (77,) 

\K h e3r h J 


where, in the last estimate, we have made use of the second estimate in Lemma 4.1 For the edge 
integrals in 9§h involving the upwind flux, we first observe that 


We then have that 


’eh 


\{w h -,n h }{ri}\ = |{w /j ;n, ! }||{77}| <^J{tj}|. 

{w h -,n h }{ri }[*]ds h < J^ h { 77 } \f^ h [x] 

L 2 (e h ) 


L 2 (e h ) 


Combining this with 


/ &*fa][*]ds h < \f^ h [n] \[$e h {x\ 


’ e h Lfe h ) 

we may bound the flux integral term of 98h as follows: 


L 2 fh) 


E / ({w/,;n/ ! }{77} + ^j77])[x]ds h < E ( + \f^M ) JieM 

\ v L 2 (e h ) V tK..\ 1 v 


e h e<S h Je h 


ei,e<S‘h 


Ll (e h )j 


L2 (e h ) 


<h 3/2 M H 2 (r) \^ L , ( J 


\ 1/2 
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where, again, we have made use of the second estimate in Lemma 4.1 to obtain the last estimate. To 
obtain estimates for (v)> we first note that 

[w h \n h \ = [P h w h -w~ l -,n h ] + [w~ l -,n h \ = [P h w h -w~ l -,n h ] + [w~ L ,P~ l n h -n~ 1 ]. 

And so, by assumption ( |3.6| l and the geometric estimate ( |2.8[ ), we have that 

\\[wh\n h \\\ L ~ ( g h) <h Pl +h 2 . 


Hence, 


E / 7Z}<IIK;»/ i ]||^(4) E \\{ri}\\me h )\\{x}\\me h ) 


e h^h &h ^ e h^h 

< {h Pl +h 2 )h 3 ^ 2 \\u.\\ H 2^h~ l ^ 2 \\x\\ L 2^ 

<(/^ +1 +/* 3 ) \\u\\ H2{r) \\ X \\ L 2 { r h y 

The first two terms of 54 (given in ( ]3. 1 0| >) can be easily shown to scale like the desired final estimate by 
applying the projection estimate given in Lemma |4~2| and making use of the second estimate in Lemma 
4.1 as before. For the last term in (3.10|>, we will require the following inverse estimate, adapted from 


Lemma 4.4 in Antonietti et a/.~]( |2014| ): 

l|Vr,zlL 2 (^)<^ 1/2 ||VrallL 2( ^)- 

Making use of this estimate, we have that 


E / ds h< E PeC'W 

e h e£h Jeh e h e£ h 

tV 2 r 


Ll ( e h) 


13 eh 1/2 {eS/r h X^h} 


L 2 (e h ) 


L 2 (dK„) 


< E MM 

K h e&k 

<he l / 2 \\u\\ H 2 {r) \\\ X \\\. 


l/2 Vr h X I 


L 2 (,K h ) 


□ 

For the first term on the right-hand side of (4.31, we require the following perturbed Galerkin or¬ 
thogonality result: 


Lemma 4.4 Let u £ H s (r), s 4 2, and u/, £ V/, denote the solutions to ( |4. 1 1 and (3.2 1 , respectively. We 
define the functional £), on V l h by 


Eh{vf,) := £?h{u ' - ) = Io + h+If 


where the three contributions Eq : Eb, and Ef meassure the variational crime from the elliptic, the hy¬ 
perbolic, and the right hand side descritization, respecively. 

Then if w/, satisfies ( |3.5| >and ( |3.6[ i, we have that 

|/ D (v ft )| + |/ / (v,)|</t 2 ||/|| i2(r) |||v A ||| , (4.5) 

\lB(y h )\ < ( h 2 + h Pl +s w h~ 1 (h 2 + h Pl ))\\u\\ H 2( r - ) \\v h \\ H , (4.6) 

for h small enough. The constant s w is defined in Theorem |4.4| 
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Before we give its full proof, we will complete that of Theorem 4.1 assuming this result. Starting 
again with the splitting (4.2 1 and (|4.3|>: 


(-«/, III < |||m / -M/il 


< 


\n\ 


[^h{u-u l h ,x)-^h{n,x) 


< 


2 + \\u\\ H 2 + [h 2 + £ 2 hJ \\f\\ L 2 + (ji£ 2 +h 2 ^J \\u\\ H 2 


Which completes the proof of our main Theorem. 

Proof of Lemma p~4| Using the definition of m/, and the fact that a solves the variational problem ( |3.1| > 
we have 


(y h )=^h(u 1 - u h ,v h ) = #/ h (u l ,v h ) - / f h v h + / fv l h - 3 ?(u,v l h ) 

Jr h Jr 

= @ h (u l ,v h ) - @(u, v ' h ) + 28 h (u l , v h ) - S 8 (u.v l h ) + fv' h - f h v h = I D + I B + If 


For If we can simply write using that f), = j) 


-/. 


If = 


J r i}-^)fA<h 2 \\f\\ L 2 {r) 


The error coming from the diffusion part was throughly studied in Dedner et al. (2013) where it was 
shown to scale satisfy the desired bound. Finally, we rewrite the error term comming from the advection 
discretization: 


h 


Y / ~ W I’ U l ' Vr h Vh + {c + Yh)u l v h dA h + Y / w h u l [n] +-[w h -,n h )u 'v h ds h 
K h e3r h jKh e h eS h Je h 

+ tf Wh {u~ l ,v h ) - Y / , -wu-V r v l h +cuv l h dA- Y / / ds 

K h e3T h jK h e h e S h Je h 


Noting that u, u 1 are continuous functions the numerical fluxes reduce to 

wijii ^ = {wh',nh}u~ l , wu = {w;n}u 

Next we use the integration by parts formula ( |3.1| > again taking into account that [ u} = [u ; ] = 0, arriving 
at: 

h= Y [ w hVh-Vr h u~ l + (c + r h +Vr h -w h )u~ l v h dA h - Y [ \[ w h\n h ]u~ l {v h } 
K h esr h jK h e h eS h Je h z 

— Y i wv h ' Vrw + cm{ dA 

K h Kh 


Note that for h small enough y/, + Wp h ■ w/ t = 0 based on the definition (|3. 3 [> and using that • u h —> 0. 
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We now lift the volume integrals on I/, to F 

h= Y f 8k l w l h v l h -Ph{I-dH)PV r u + cuv l h dA- Y f, \[wh\n h ]u~ l {v h } 

K h£^h K h e h 6 <4 4 

- E .wv l h -Vru + cuv l h dA 

= E [ ,(8 h -l)(w l h v l h -P h (I-dH)PVrudA 

K h e.% JK h 

+ Y [ I ( w h - w ) v h ■ p h ( 7 - dH)PV r u dA 

K h e3r h jK h 

+ E /, ( wv l • ( p * ( 7 - - 7 ) Vr« dA 

K h e,9 h jK h 

~ Y f, ^[ w hi n h]u '{Vh} = + 7 b + 7 b + 7 b 

«*e<4 ■ /e * 7 

We will bound each of these terms using results from Lemma ( |2.1[ i and Assumptions ( |3.5[ )and ( |3.6[ > : 

4 < h 2 H\\ L ~ {n \\Ph(I-dH)P\\ L ~ {r) Y /. bX2 { Kl)\\M\ L 2 ( Ki) £ h 2 \\v l h \\ H \\u\\ H 2 {r) 

K h e3r h JK h 


ll< Y , \\vX^Ki)\\^ l h-^)-Ph{I-dH)P^ r u\\ <h^\\v l h \\ H \\u\\ H2{r) 

K h e3r h JK i, 

For the next estimate we use that w is tangent to F so that w = Pw and that P is symmetric: 

4 = E [, (v l h w{PP h {I-dH)P-P)V r udA 

K h e3r h jK h 

= E [ l vl h w -( pp h p -P)^ru+dw-P h HPV r udA<li 2 \\v , h \\ H \\u\\ H 2 (r) 

K h e3r h jK h 

It remains to bound the term on the skeleton of the grid. Note first that 

Ki;»/<] = (« - W~ l ) -n+) + (O^ -w~ 1 )- n h ) + W -'- (n+ + n h )< h Pl +h 2 

using our Assumptions on w/, and scalling results from Lemma ( |2. 1 [ i. This leads to the following esti¬ 
mate 

4 < (h pi +h 2 ) Y h~ l \\L^e h) \\n\\ LHeh) < 0 h p1 +/. 2 )/.- 1 ||«- , || tf(rA) ||v*|| tf(Ii) 

eh&^h 

This completes the proof. □ 


5. Construction of discrete velocity field 

We will now attempt to justify the assumptions we have made on w/, by constructing a discrete velocity 
field which satisfies assumptions (3.5) and ( |3.6| l. We will first discuss why we do not simply take 
Wh = w~ l and then describe two alternative approaches, 
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5.1 Downward lift of velocity field 

Consider the simplest choice for a discrete velocity field given by w\, := vv 1 . Note that due to the 
definition of our bilinear form we can always write P\ x w\ x instead of w/,, so using the projection onto the 
tangent planes of the triangles is not required in the definition of the discrete velocity field. 

In general one can neither expect this choice to lead to a divergence free field on each triangle, nor 
that the normal jumps across edges will vanish, i.e., [w^nf = +nf)f=0. We added a number of 

terms to take into account that the velocity field w/, is not divergence free and has non continuous normal 
components over the element edges. Simply defining the discrete bilinear form by taking vv/, = vv 1 in 
f'/ih, dropping the extra c to Wh term and taking jh = 0 does not lead to a unconditionally stable scheme. In 
fact the matrix resulting from such a scheme will not be positive-definite independently of h. To see this 
consider the bilinear form 

Y / -w~'uh ■ Vr A v ft + U h v h dAh + Y / w ~ lu h[vh\ dsh 

K h e3T h jK h e h eS h Je h 

Integrating by parts and using Lemma 0 as in the proof of Lemma |3.2| choosing «/, = 0 for every 
Kh £ ,% except for two elements and Kf for which e/, = K h fj K h . Furthermore, we can assume 
without loss of generality that nf = (—1,0,0), rif = (cos^), sin(<gr),0) with q £ (0,2tr). Note that 
unless q = 0,27f, we have that nfi —nf. The velocity w~ l at e/, is assumed to be (—1,0,0), so that 
w _/ • = 1 > 0 and w^ 1 ■ nf = — cos (q) < 0. Finally, we assume that u^ = uf = 1 so that [«/,] = 0 on 

e/,. With these conditions, the stability of ( |3.2[ > boils down to showing that 

Y L (w~ l ■ n Kj \ uj, ds h ^ 0 

z K h e.% Jr)K h K ' 

which, from Lemma [3TT| and the above conditions, is equivalent to showing that 

f [w~ l Uh\n h ] ds h ^ 0 . 

Notice that the numerical flux does not appear given that it is scaled with [uf\ = 0, and thus cannot 
influence the sign of the above quantity. Expanding the expression, we have that 

If i r 1 

-- / [w~'u 2 h ;n h ] ds h = -- / w~ l u 2 h + -n£ + w~'u 2 ~ n h ds h = -\e h \(cos(q) - 1 ) < 0 . 

^ J eh ^ ^ 


Hence, in general, whenever f —nj, /i-indepcndent positive-definiteness of the matrix resulting from 
the scheme may not hold, regardless of the choice of the modified upwind flux. 

So in fact both Jh and r S Wh are important terms to make the scheme positive definite independent of 
h. Since the ^Wh will not vanish, our error estimate indicate a suboptimal convergence rate (of course p\ 
is arbitrarily large), although as already pointed out, this is not confirmed by our numerical experiments. 
A more severe problem is, that evaluating //, requires the computation of V^ • vv 1 , which requires 
derivatives of the lifting operator and thus more information about the surface F then we wish to have 
in our numerical scheme. We therefore will not use this choice in our numerical experiments, but tests 
indicate that this choice is comparable to the other choices described in the following. 
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5.2 Lagrange interpolation 

We can use a Lagrange interpolation of w~ l on I/, to define w/,. This is easy to implement and evaluating 
both y/, and ^w h is not problematic. According to our error analysis the approximation order of linear 
finite elements is sufficient (p\ = 2 ) but of course % Vh will in general not vanish and thus our error 
estimate will be suboptimal. In our numerical experiments we still observed an optimal rate while 
a piecewise constant interpolation does not lead to optimal results. We omitted the detail of these 
experiments in this paper. 

5.3 Surface Raviart-Thomas interpolant 

Our next choice avoids the problem of suboptimality by constructing a velocity field with = 0. 

To this end we make use of a Raviart-Thomas-type interpolant of w f which we will refer to as the 
surface Raviart-Thomas interpolant. 

Let F Kh denote the mapping from the reference element K to Kj v Then we have that V F Kh = (eo,ej ) F 
R 3x2 where eo and e\ are two edges of Kj, intersecting at the vertex xo- We first define the local spaces 

V l RT {K h ) := {s h (x) := VF Kh (F^(x))p (/£(*)) , p £ pP«(tf)] 2 } . 

We next define the local Raviart-Thomas space of order q on K/, to be given by 

RT q (K h ) ■= jw A (x) := s h (x) + (x-x 0 )t (^(x)) , s h £ P % r (K h ), t £ . 

It is clear from the definition of RT q (Kh) that any function w/, £ RT q (Kf) for every K /, £ is tangential 
to i~7z- Using the convention that the conormal to eh C dKh is nt, the local degrees of freedom of 
Wh £ RT CI (K/,) are given by 


' e h 


w h • rifi P q dsh G F q (e h ),e h cdK/,, 


[ W,, ■ P q - 1 ds h Vp 9 _ tePj/ (K h ). 

J Kh 


(5.1) 

(5.2) 


We then define, for w 1 £ \W 1 "(W/,)] 3 , the local surface Raviart-Thomas interpolant of order q to be 
n q K w-‘ £ RT q (Kh) satisfying 






-l 




-l 


■ Ilf p q ds h = / w 1 ■ rif h p q ds h Wp q £ P q (eh),eh C dK h , 

Je h 

■Pq-\ ds h = / w~ l -p q -i ds h Vp f/ _i £ P It 1 (Kh). 

JK h 


(5.3) 

(5.4) 


Here, the “average” conormals n e f are given by n e f := ±-j—v—— 

\2\ n h~ n h'\ 

Remark 5.1 Notice that this definition differs from that of the local classical Raviart-Thomas inter¬ 
polant in the way we have defined the right-hand side of (|5.3[). We have to use what we call the “average” 


conormals nf[ instead of the standard conormals nf because they satisfy «+ = —n e . From here on, 
we will refer to the local classical Raviart-Thomas interpolant by 17^ vv 1 . 


+h 
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Lemma 5.1 Let 771% w 1 be the local surface Raviart-Thomas interpolant of w 1 £ [W 1, “(.K'a)] 3 , defined 
as in (|5.3|)-(|5T4|>. Then we have that 


n K h w ' n h = ~ n k h W ' n h 


for each eh £ S'h- 


Proof By (|5.3|»—(|5.4)» and using the fact that nj = —n eh , we have that 


/ n K h w l - n hPq ds h = / W l -n+ h p q ds h = — / w 1 ■ n eh p q ds h = - / 17% w 1 ■ n h p q ds h . 
Je h Je h Je h J e h 


It follows that 


for every p q £ f* 1 (e/,). By Proposition 3.2 in 
which gives us the pointwise equality 


[ ( n l, w '■ n h + n l„ w l -”h) Pq ds h = 0 

Je h v ' 

11991 f, we have that 77% i w _/ • n% £ P 9 (e/,), 


Fortin & Brezzi 


n K h w ’ n t = - n K h w ■ n h 


as required. □ 

Lemma 5.2 Let 77% i w~' be the local surface Raviart-Thomas interpolant of w~ 1 £ [W I,00 (X*)] 3 defined 
as in (5.3 M 5.4 1 and let 17% w be its local classical Raviart-Thomas interpolant. We then have that 


\\K w ~ l -n q Kh 


W l \\L°°{K h UdK h ) ^h 2 


for each Kj, £ 3f t . 
Proof. Denote by 


dK n dK 

* the set of local degrees of freedom given by (5.1 1 and {(p. h } i= % 

j K h \ nK h 


the 


associated (vector-valued) basis functions. Similarly, we denote by {/V ( - *} . \ the set of local degrees of 

K riff 

freedom given by (5.2 1 and ( <p i h } i= \ the associated (vector-valued) basis functions. The local degrees 


of freedom for the local standard Raviart-Thomas interpolant { N, }. 
We then have that 


n 3K,, 


and {Ni} i= \ are defined similarly. 


n 3K u 


n Ku 


n q K w~ l (x) = £N? Kh (w- l )(p? Kh (x) + £N?>'(w- l )(p? h (x), 
1=1 i=i 


and similarly for 77^ w l . Then by noting that Nf h (w l )=Nf h (w l ) and making use of (5.1! and (5.31 
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we have that 


II n K h w l ~ n K h W l \\L°°(K h UdK h ) 


n BK h 

1=1 v ' 


L°°(K h UdK h ) 


f max 


= max 
1 <i< n 3K h 


< max 
l*H<n SKh 


N? Kh (w- l )-N? Kh {w~ l ) 


dK/j 

N- Kh (■ w ~ l ) - ivf Kh {w~ l ) E 

i= 1 

[ w ~ l ■ n t h & ds h - [ W~ l ■ tl+ & ds h = , max / w_/ ■ ( P ~ V* - P ~‘ n h) & ds h 

Je h Je h l<i < n 3K h Je h ' ' 


^ ll«- p »J"JlL“-(4) + ll n “ P ”/IIU°“(4) - /?2 


2.1 


□ 


basis functions of P 9 (<?/,). The last estimate follows from Lemma iz. i lj 

^ <u -'- the case of the local surface Raviart- 


where {<!;,■} denote the basis functions of P q (eh). The last estimati 
” " " ving theorem will help justify assumption (3.5 ifor 

uaI finf nf 7Prn nrH<=»r (n — 


The following tneorem win netp jusn: 

Thomas interpolant of zero order (q = 0). 

Theorem 5.1 Let w“ z G [W l =“(A'/ ! )] 3 and il£ w 1 be its local classical Raviart-Thomas interpolant of 
zero order defined only through condition (5.31 (with replaced by nff). We then have that 

HP/jVV ’ ~n^ h W l \\ L ~ {Kh) <h\\Vr h -W ( || L°°(K h )i 


Vn- 


( P AW" 




'Kt, 


for each Kh G fTh. 


Proof. The proof of the first estimate follows similar lines as that of Theorem 6.3 in Acosta et al. 
( [2011) >. The second estimate follows similar lines as that of Theorem 1.114 in jErn| ( |20Q4[ ). □ 

The first estimate of Theorem |5 .1 | together with Lemm a |5,2| gu arantees that the local surface Raviart- 
Thomas interpolant also satisfies Theorem |5.1| As such, assumption (3.5 1 holds when choosing w/ ; to 
be the local surface Raviart-Thomas interpolant of zero order. 

We finally show that assumption (3.6( holds for local surface Raviart-Thomas interpolants of zero 
order. 


Lemma 5.3 Let w 1 G [W 1 “(A - /,)] 3 , Kh G , and 172 vv 1 be its local surface Raviart-Thomas inter¬ 
polant of zero order defined only through condition (|5.3|). We then have that 




< h. 

L-(K h ) 


Proof. We have that 


V^.< vv"' 


L°°(K h ) 


€ 


V L ■ { n K h w l - p hW~ 


|Vr A - p *w 


L~(K h ) 


Making use of Lemma 3.2 in Olshanskii et al. ( 2013] ), we have that the second term scales like h. For 
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the first term, we have that 


Vr*-( 


n° Ki w- l -p„w~ l 


L°°(K h ) 


< 


Vr*-( 

Vr*-(. 




17^ w z —P/,w ( 




£"(**) 


The second term in the above scales appropriately by the second estimate of Theorem 5.1 For the first 
term we proceed as in the proof of Lemma [5^2] to get that 




■ (< w-'-Mvw- 




^ max 

< max 
KKn dKh 


N- K \w~ l )-N- Kh (w~ l ) £ V L • % 


L°°(K h ) 
dK h 


dKu 


N- Kh {w- l )-N™ h {w~ l ) 


■dK h (...-b 


1 <:i<n 3Kh 

hr 1 < /j 


as required. 


□ 


6. Numerical tests 

For the test problems discussed below, we will focus on a surface IP discretisation of the diffusion term 
and call the resulting approximation the surface IP/UP approximation. Furthermore, the discrete velocity 
field w/, is chosen to be the zero order surface Raviart-Thomas interpolant of vv / i.e. w/,| K h = JTj^ vv ; . 
We will also briefly discuss the case when we choose vv;, = vv ; in the numerics. 


6.1 Test problem on torus 

Our first test problem, considered in Olshanskii et al. ( [2013] ), involves solving ( |3. 1 ] ) on the torus 

r= |(xi,x 2 ,x 3 ) | + X 3 = J^ 

with velocity field 

1 T 
w(x) = (—x 2 ,xi,0) . 

Note that the velocity field w is tangential to the torus and divergence-free. We set e = 10 6 and 
construct the right-hand side / such that the solution u of is given by 

, N *1*2 . f X3 

mix) =-arctan 

K 

Note that u has a sharp internal layer as shown in Figure [l] 

Figure [I] shows the exact solution and both the unstabilised surface FEM approximation and the 
surface IP/UP approximation of Notice how, as in the planar case, the unstabilised surface FEM 
approximation exhibits global spurious oscillations whilst the surface IP/UP approximation is com¬ 
pletely free of such oscillations. We obtain similar results for the case when we choose vv/, = vv 1 in the 
surface IP/UP method, although L°° errors tend to be slightly larger for such a choice. 
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Fig. 1: Exact solution of (3.1 1 (top) and pointwise errors for respectively the (unstabilised) surface 
FEM approximation (bottom left) and the surface IP/UP approximation (bottom right) on the torus 
(1410 elements). 
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6.2 Test problem on sphere 

Next, we consider ( |3.1[ ) on the unit sphere 

r = {rgR 3 : |x| = 1} 


with velocity field 

w(x) = ^ X2 1 1 —*3,0 

Notice again that w is tangential to the sphere and divergence-free. We set e = 10 ( ' and construct the 
right-hand side / such that the solution u of is given by the expression given in the previous test 
problem. Tables [I] and [5] show the L 2 and DG norm errors/EOCs outside the sharp internal layer, given 
by D = {x £ r : |x 3 1 > 0.3}, for respectively the (unstabilised) surface FEM approximation and the 
surface IP/UP approximation. 



Elements 

h 

L2(D)-error 

Z/ 2 (.D)-eoc 

DG(D)-e rror 

DG(D)-e oc 

632 

0.2239 

0.04462 


0.865 


2528 

0.1121 

0.01736 

1.36 

0.652 

0.40 

10112 

0.0561 

0.00936 

0.89 

0.727 

-0.16 

40448 

0.0280 

0.00604 

0.63 

0.934 

-0.36 

161792 

0.0140 

0.00356 

0.76 

1.095 

-0.23 

647168 

0.0070 

0.00169 

1.07 

1.038 

0.08 


Table 1: Errors and convergence orders for the (unstabilised) surface FEM approximation of (3.11 on 
the subdomain D of the unit sphere. 


Elements 

h 

L2(7J)-eiTor 

Li (D)-eoc 

DG(D )-error 

DG(D)-e oc 

632 

0.2239 

0.0073256 


0.15932 


2528 

0.1121 

0.0021745 

1.75 

0.08892 

0.84 

10112 

0.0561 

0.0006499 

1.75 

0.05015 

0.83 

40448 

0.0280 

0.0001917 

1.76 

0.02820 

0.83 

161792 

0.0140 

5.399e-05 

1.83 

0.01537 

0.88 

647168 

0.0070 

1.394e-05 

1.95 

0.00778 

0.98 


Table 2: Errors and convergence orders for the IP/UP approximation of (3.1 1 on the subdomain D of the 
unit sphere. 


The results clearly indicate that the surface IP/UP method performs better than the unstabilised 
surface FEM. The results for the surface IP/UP method indicate a 0(h 2 ) convergence in the Lr(D)- 
norm and 0(h) in the DG(D)-norm. The unstabilised surface FEM, on the other hand, shows a much 
more erratic behaviour and does not attain its asymptotic convergence rates within our computational 
domain. 

Table [ 3 ] show the relevant errors when using w/, = w 1 in the surface IP/UP approximation. The 
errors appear to be smaller by a factor of about 0.5 compared to those shown in Table [2] for which 
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Elements 

h 

L2(7l)-error 

Li(Z))-eoc 

DG[D) -error 

DG[D)-t oc 

632 

0.2239 

0.0040846 


0.11275 


2528 

0.1121 

0.0010464 

1.96 

0.05707 

0.98 

10112 

0.0561 

0.0002654 

1.98 

0.02867 

0.99 

40448 

0.0280 

6.679e-05 

1.99 

0.01437 

1.00 

161792 

0.0140 

1.670e-05 

2.00 

0.00718 

1.00 

647168 

0.0070 

4.161e-06 

2.00 

0.00359 

1.00 


Table 3: Errors and convergence orders for the IP/UP approximation of (3.1 1 with Wh = w 1 on the 
subdomain D of the unit sphere. 


we chose Wh\K h = 1 . This can be explained by the fact that triangulations for simple surfaces 

such as the unit sphere can be constructed to be very “smooth” (in the sense that the relation nt = 
—nj t practically holds for each e/, £ Sj t ) and that the zero order Raviart-Thomas approximation error is 
relatively large. 

7. Conclusions 

It is well known that the DG method is especially well suited for stabilizing transport terms in PDE 
models. This type of problem has not yet been studied using discrete surface finite-elements. In this 
paper we proved a-priori error estimate for the DG method for the stationary linear hyperbolic problems 
and elliptic problems with possibly dominate advection term. This extends previous results for the 
laplace equations on surfaces. The theory and the numerical experiments show that special care has 
to be taken when projecting the given continuous velocity field unto the discrete surface. The main 
two problems identified are the non zero divergence and the jump in the normal component. Different 
suggestions were discussed on how to handle these problems either by extending the bilinear form to 
take the jump in the velocity field into account or by the use of a special projection based on Raviart- 
Thomas like interpolation operators. Numerical experiments demonstrate the accuracy of the resulting 
method. 
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